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ABSTRACT 

This  thesis  examines  acoustic  transient  discrimination  and  Time  Difference  Of 
Arrival  (TDOA)  estimation  for  the  purposes  of  estimating  the  position  of  a  submarine  in 
a  sonabuoy  field.  Transient  discrimination,  for  this  thesis,  is  the  process  of  telling 
different  transients  apart.  Two  algorithms  are  evaluated.  One  method  is  based  on  higher 
order  statistics  while  the  other  is  based  on  signal  subspace  techniques.  Extensive 
simulations  using  synthetic  transients  were  conducted  to  establish  the  performance  of 
each  algorithm  in  terms  of  discrimination  and  TDOA  estimation.  It  was  found  that  the 
bispectral  algorithm  gave  better  TDOA  estimation  at  low  SNRs  while  the  subspace 
algorithm  gave  better  TDOA  estimation  at  high  SNRs.  For  discrimination,  it  was  found 
that  the  subspace  algorithm  gave  consistant  false  alarm  rates  at  all  SNRs  while  the  false 
alarm  rate  for  the  bispectral  algorithm  grew  with  increasing  SNR. 
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EXECUTIVE  SUMMARY 

In  this  thesis  we  have  developed  and  compared  two  algorithms,  namely  the 
bispectrum  and  subspace  linear  phase  detectors.  These  algorithms  were  developed  for 
the  purposes  of  transient  discrimination  and  Time-Difference-Of-Arrival  (TDOA) 
estimation.  They  are  to  be  used  as  part  of  a  transient  tool  suite  to  aid  in  the  estimation  of 
a  submarine's  position.  Two  performance  measures  were  used  to  evaluate  the 
algorithms,  namely  the  probability  of  correct  TDOA  (Pj)  and  the  probability  of  correct 
discrimination  (PdD- 

In  general,  it  can  be  said  that  for  TDOA,  the  bispectral  linear  phase  detector  gave 
better  results  at  low  SNRs  while  the  subspace  linear  phase  detector  worked  better  at  the 
higher  SNRs. 

For  discrimination,  it  was  found  that  the  bispectral  discriminator  gave  higher  Pdi 
than  the  subspace  discriminator.  However,  the  probability  of  false  discrimination  (Pfa)  of 
the  bispectral  discriminator  increased  at  higher  SNRs  while  the  subspace  discriminator 
gave  a  constant  P&  for  all  SNRs  evaluated.  For  discrimination,  the  advantage  of  having  a 
constant  Pfa  is  desirable.  Therefore  the  subspace  discriminator  is  the  best  option  even 
although  it  produced  lower  Pdj  than  the  bispectral  discriminator.  It  was  also  found  that 
there  are  design  trade-offs  between  processing  speed  and  performance  that  need  to  be 
made.  For  the  bispectral  linear  phase  detector,  this  trade-off  is  in  terms  of  threshold  gain; 
for  the  subspace  linear  phase  detector,  this  trade-off  is  in  terms  of  correlation  matrix  size. 
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I.    INTRODUCTION 

The  work  of  this  thesis  forms  part  of  the  ongoing  effort  to  automate  the  detection, 
discrimination  and  Time-Difference-Of- Arrival  (TDOA)  estimation  of  transient  signals. 
At  present,  transient  signals  are  located  and  processed  manually,  making  it  a  very  time 
consuming  procedure.  It  is  therefore  more  desirable  to  have  an  automated  system  that 
can  deliver  the  same  or  better  results  than  the  present  manual  approach. 

For  this  thesis,  it  is  assumed  that  the  detection  of  the  transient  has  already  taken 
place;  therefore,  the  objective  of  the  thesis  is  to  look  at  discrimination  and  TDOA 
estimation  only.  Two  algorithms,  namely  the  bispectral  linear  phase  detector  and  the 
subspace  linear  phase  detector,  are  developed  for  this  purpose. 

A.         THE  OPERATIONAL  MISSION 

A  typical  mission  considered  for  this  thesis  starts  off  with  an  anti-submarine 
warfare  (ASW)  aircraft  laying  a  sonabuoy  field  in  the  vicinity  or  on  the  predicted  path  of 
a  submarine.  Usually  each  sonabuoy  in  the  field  consists  of  calibrated  omni-directional 
passive  acoustic  sensors.  The  ASW  aircraft  collects  acoustic  data  on  a  target  submarine 
as  it  passes  through  the  sonabuoy  field. 

Figure  1  shows  a  typical  V-shaped  sonabuoy  field  that  can  be  used  in  this  type  of 
mission.  When  this  shape  is  used,  the  target  must  pass  through  the  apex  of  the  V  of  the 
field  (as  close  as  possible)  in  order  to  obtain  the  most  accurate  results.  The  data  received 
at  the  sonabuoys  is  transmitted  to  the  aircraft,  which  in  turn  records  the  data  on  tape  for 
mission  post  processing. 


Figure  1.  Typical  Sonabuoy  Field. 

The  data  collected  on  the  tapes  is  processed  at  an  onshore  location  in  order  to  estimate  the 
acoustic  Sound  Pressure  Levels  (SPL)  of  tonals,  broadband  signals  and  transients  that  are 
emitted  from  the  submarine. 

To  obtain  precise  SPL  of  a  submarine,  an  accurate  track  is  required  so  that  the 
pressure  levels  received  at  the  sonabuoys  can  be  projected  back  to  estimate  the  pressure 
levels  at  the  target.  For  best  results,  the  position  of  each  sonabuoy  must  be  known,  the 
target  must  be  on  a  course  that  passes  through  the  middle  of  the  field,  and  an  accurate 
estimation  of  TDOA  of  acoustic  signals  must  be  made. 


B. 


DATA  PROCESSING 


The  data  collected  from  the  sonabuoys  and  stored  on  the  tape  is  processed  by  first 
generating  a  baseband  sonogram.  A  sonogram  is  a  frequency  versus  time  plot  of  the 
acoustic  pressure  levels  at  a  sonabuoy.  From  these  sonograms  an  analyst  can  identify  and 
highlight  contact  events  of  interest,  which  may  include  signals,  such  as  narrowband 
tonals  and  transients.  The  events  identified  are  then  processed  to  obtain  an  estimate  of 
the  target's  track.  Once  the  best  estimate  of  the  target  track  has  been  made,  the  analyst 
can  conduct  SPL  analysis  on  all  types  of  previously  identified  contact  signals.  The  final 
desired  results  are  the  SPL  signature  characteristics  for  a  given  target  with  respect  to 
aspect  angle. 


At  present,  transient  analysis  is  a  manually  intensive  operation.  Typically,  the 
analyst  marks  transient  events  of  interest  from  the  sonogram  display  of  a  sonabuoy  and 
then  searches  the  sonograms  of  the  other  sonabuoys  for  the  same  transient.  This  can  be  a 
very  time  consuming  process  since  in  a  single  mission  there  can  be  hundreds  of  transients 
within  a  single  sonogram.  These  must  be  matched  up  to  those  on  the  sonograms  from  the 
other  sonabuoys  in  the  sonabuoy  field. 

C.  THESIS  GOAL 

The  goal  of  this  thesis  is  to  develop  and  evaluate  two  algorithms  that  can  be  used 
for  transient  discrimination  and  Time-Difference-Of-Arrival  (TDOA)  estimation.  The 
two  algorithms  are  the  bispectral  linear  phase  detector  and  the  subspace  linear  phase 
detector. 

D.  THESIS  OUTLINE 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  II  discusses 
transients  and  transient  processing  and  introduces  the  problems  of  interest,  namely 
TDOA  estimation  and  transient  discrimination.  Chapter  HI  presents  the  relevant  signal 
models  and  the  analysis  of  the  measured  signals  in  terms  of  second-  and  third-order 
moments.  This  chapter  details  the  formulation  of  the  bispectral  linear  phase  detector  and 
discriminator.  Chapter  IV  discusses  signal  subspace  techniques  and  their  application  to 
the  problem  of  TDOA  estimation  and  transient  discrimination.  The  outcome  of  this 
chapter  is  the  formulation  of  the  subspace  TDOA  estimator  and  discriminator.  Chapter  V 
examines  the  results  achieved  by  both  the  bispectral  and  subspace  TDOA  estimators  and 
discriminators.  Extensive  simulations  are  conducted  using  synthetic  transients  to  produce 
performance  curves  and  to  study  the  utility  of  each  technique  in  discrimination  and 
TDOA  of  transients.  Chapter  VI  provides  conclusions  and  recommendations  for  future 
work  based  on  the  results  presented.  Appendix  A  and  Appendix  B  show  the 
discrimination  results  for  the  different  combinations  of  signals  not  shown  in  Chapter  V. 
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II.    TRANSIENT  PROCESSING 

A.         TRANSIENT  SIGNALS 

Transients  of  interest  here  are  short  duration  wideband  acoustic  signals. 
Accordingly  transients  can  be  of  varying  shapes  and  lengths,  lasting  anywhere  from  a  few 
microseconds  to  a  couple  of  seconds.  Typical  examples  of  these  signals  are  a  wrench 
falling  on  a  metal  deck  or  the  sounds  of  a  biological  life  form,  such  as  a  whale.  Due  to 
the  diversity  of  these  signals,  there  is  limited  a  priori  information  that  can  be  used  to  aid 
in  their  detection.  When  transients  are  imbedded  in  noise,  their  detection  becomes  very 
difficult  because  the  energy  in  the  noise  frequently  dominates  the  energy  in  the  transient 
over  the  interval  of  interest.  Nevertheless,  in  conjunction  with  other  measurements  or  by 
themselves,  transients  can  be  used  to  estimate  the  track  of  the  target  submarine  by 
estimating  the  TDOA  of  a  transient  between  two  sonabuoys. 

Tracking  of  a  submarine  using  transients  requires  transient  detection, 
discrimination  and  TDOA  estimation.  Prior  to  discrimination,  a  transient  must  be 
detected  in  a  noisy  background  on  the  sonagram  of  a  single  sonabuoy.  This  thesis  does 
not  address  the  transient  detection  problem  and  therefore  assumes  that  a  transient  has 
been  detected  on  the  sonogram  of  a  sonabuoy. 

After  the  transient  has  been  detected,  the  same  transient  must  be  found  on  the 
sonograms  of  the  other  sonabuoys  in  the  sonabuoy  field.  This  is  difficult  since  the 
detected  transient  may  not  exist  on  all  sonabuoys  within  the  field.  That  is,  the  sonograms 
of  the  other  sonabuoys  could  contain  only  noise  or  even  a  different  transient  in  the 
window  of  interest.  The  first  case  (i.e.,  where  the  sonogram  of  the  other  sonabuoys 
contain  only  noise)  is  a  common  scenario.  This  happens  because  transients  can  be  very 
localized  and  thus  do  not  appear  on  the  sonograms  of  each  sonabuoy.  The  second  case 
(i.e.,  where  a  different  transient  exists  on  the  sonogram  of  the  second  sonabuoy)  is  also 
common  due  to  the  complex  environment  of  the  mission,  where  multiple  transients  from 
different  sources  can  exist  on  the  data  records  of  the  other  sonabuoys  in  the  field.  It  is 
therefore  important  to  discriminate  between  different  transients  and  to  know  if  there  is  no 
transient  present.    If  the  transients  arriving  at  two  or  more  different  sonabuoys  are  the 


same,  then  TDOA  estimates  can  be  used  for  target  localization.   In  the  following  we  first 
discuss  TDOA  estimation  and  then  transient  discrimination. 

B.         TIME  DIFFERENCE  OF  ARRIVAL  (TDOA)  ESTIMATION 


TDOA  estimation  is  by  no  means  simple  and  therefore  many  authors  have  written 
about  it  [Ref  1].  It  is,  however,  the  primary  means  of  determining  range  to  the  target  in 
passive  detection,  since  TDOA  information  from  multiple  sonabuoys  and  the  geometry  of 
the  sonabuoy  field  can  be  used  to  determine  the  position  of  the  submarine  at  any 
particular  instance  in  time.  For  any  given  value  of  TDOA  between  two  sonabuoys,  the 
locus  of  possible  target  positions  is  a  hyperbolic  curve.  If  more  than  one  curve  can  be 
drawn,  i.e.,  if  TDOA  values  between  three  or  more  sonabuoys  can  be  calculated,  the 
intersection  of  the  curves  determines  the  position  of  the  target. 

Consider  a  simple  two-buoy  model  as  shown  in  Figure  2.  The  sonabuoys  are 
symmetrically  positioned  about  the  origin  on  the  x-axis  as  shown  in  Figure  2,  with  the 
target  located  at  T. 


Buoy  2  (-x0,0) 


Buoy  1  (xo,0) 


Figure  2.  Two  Symmetrically  Positioned  Sonabuoys. 


For  this  geometry  the  time  delay  r  between  sonabuoy  1  and  sonabuoy  2  is  given 


by 


T=t.-t0  = 


A,     Dh 


(2.1) 


(r  •  cf  =  (-y](x0-x)2+y2  -  ^](x0+x)2  +  v2  )2 


where  c  is  the  speed  of  sound,  D,    and  Dt    are  the  distances  between  the  target  and 

sonabuoy  1  and  sonabuoy  2,  respectively.  After  some  algebraic  manipulation1  of  Eq.  2.1, 
the  following  expression  is  obtained  [Ref  1]: 

x 


a 


where 


2  b- 


T.C 


(2.2) 


b  = 


2xr 


—  a 


Equation  2.2  describes  a  hyperbola,  and  a  set  of  hyperbolas  can  be  drawn  for  different 
TDOA's  between  sonabuoy  1  and  2  as  shown  in  Figure  3. 


Figure  3.  Hyperbolic  curves  for  TDOA  between  two  symmetrically  positioned 

sonabuoys.  Sonabuoy  positions  are  indicated  by  the  dots. 


C. 


TRANSIENT  DISCRIMINATION 


The  objective  of  discrimination  between  transients  is  to  be  able  to  tell  different 
transients  apart.  This  is  different  from  transient  classification  or  transient 
characterization,  which  attempts  to  identify  the  source  of  the  transient  and  seeks  to  group 


similar  transients  together. 


To  obtain  a  hyperbolic  expression,  Eq  2.1  must  be  squared,  and  all  the  cross  terms  must  be  left  on  the 
right  hand  side  and  all  other  terms  on  the  left  hand  side.  Both  sides  must  now  be  squared  again  with  like 
terms  collected. 


Discrimination  is  difficult  since  it  relies  on  detecting  differences  between  the  two 
received  signals.  These  differences  can  be  either  in  magnitude  or  phase  or  both.  The 
problem  is  further  complicated  by  the  presence  of  noise  and  the  fact  that  (due  to 
differences  in  propagation  path  characteristics)  the  same  transient  arriving  at  two 
different  locations  may  not  look  and  sound  the  same.  Achieving  good  discrimination  at 
low  SNR  values  is  a  challenging  task. 

Figure  4  shows  some  typical  cases  in  which  transient  discrimination  has  to  take 
place.  The  cases  shown  in  this  figure  are  a  localized  transient  such  as  an  expanding 
bubble  on  a  sonabuoy,  a  directional  transient  that  is  emitted  from  the  hull  of  the 
submarine  and  is  only  received  at  some  of  the  sonabuoys  in  the  sonabuoy  field,  and  an 
omnidirectional  transient  that  is  received  at  all  sonabuoys  in  the  field. 


O 


Local  Transient 


o 


o 


o 


o 


.  Omni-Directional 
Transient 


o 


o 


o 


o 


"  Directional  Transient 


Figure  4.  Discrimination  Cases. 

Considering  these  cases,  the  TDOA  cannot  be  estimated  for  all  combinations  of 
two  sonabuoys  in  the  field.  Further,  if  a  transient  arriving  at  one  sonabuoy  is  mistaken  to 
be  the  same  as  the  transient  arriving  at  another  sonabuoy,  the  resulting  TDOA  estimate 
will  produce  an  erroneous  target  position.  This  in  turn  may  lead  to  errors  in  tracking  and 
ultimately  possible  erroneous  SPL  calculations. 


III.    MOMENT-BASED  SIGNAL  PROCESSING 

Most  of  the  traditional  work  on  transient  processing  and  TDOA  estimation  is 
based  on  the  use  of  second  order  statistics  and  classical  (i.e.,  Fourier-based)  methods  of 
spectrum  estimation  [Ref  2],  [Ref  3],  [Ref  4],  [Ref  5],  [Ref  6].  Recently,  new  work  has 
appeared  using  techniques  involving  higher  order  moments  of  signals  and  higher  order 
spectra,  such  as  the  bispectrum  [Ref  7]. 

The  motivation  for  using  the  bispectrum  in  transient  discrimination  and  TDOA 
estimation  is  twofold.  First,  the  higher  order  spectra  suppress  Gaussian  noise  processes 
of  unknown  spectral  characteristics.  Secondly,  these  spectra,  unlike  the  usual  power 
density  spectrum,  preserve  phase  information  [Ref  7].  In  the  last  few  years  there  has 
been  a  considerable  amount  of  new  research  done  in  using  the  bispectrum  and 
trispectrum  for  transient  detection,  time  delay  estimation  and  classification  of  signals 
[Ref  7:p  313],  [Ref  8],  [Ref  9],  [Ref  10],  [Ref  11],  [Ref  12]. 

This  chapter  defines  the  signal  model  used  in  this  thesis  and  the  analysis  of  these 
signals  based  on  second  and  third  moments.  Finally,  a  set  of  algorithms  for 
discrimination  and  TDOA  estimation  using  third  order  moments  is  discussed. 

A.         SIGNAL  MODEL 

Before  any  further  analysis  can  be  done  it  is  necessary  to  develop  a  signal  model, 
which  can  be  used  for  the  analysis  of  transients  arriving  at  two  sonabuoys.  As  mentioned 
previously,  there  are  three  cases  to  be  considered.  In  the  first  case,  the  transient  arriving 
at  the  second  sonabuoy  is  the  same  as  the  transient  arriving  at  the  first  sonabuoy.  In  the 
second  case,  the  two  arriving  transient  signals  are  different  while  in  the  third  case  only 
noise  is  present  at  the  second  sonabuoy.  It  is  also  assumed  that  the  sonabuoys  used  are 
omni-directional  and  that  their  separation  in  the  sonabuoy  field  is  large  enough  so  that 
noise  at  the  first  sonabuoy  is  uncorrelated  with  noise  at  the  second  sonabuoy  in  both 
space  and  time. 

The  following  simple  model  can  be  used  to  represent  a  single  transient  arriving  at 
two  different  sonabuoys: 


Xi(k)  =  s(k)+TJl(k) 

x2(k)  =  As{k-L)+rj2{k\ 

where  X;(k)  is  the  noise-embedded  signal  arriving  at  sonabuoy  i,  s{k)  is  the  transient 

signal  itself,  and  r/^k)  is  white  gaussian  noise.    Note  that  the  signal  at  sonabuoy  2  is 

subject  to  a  relative  attenuation  A  and  delay  L  with  respect  to  the  signal  at  sonabuoy  1. 

The  frequency  domain  expression  for  these  signals  is 

XAo))=S(o))+N.(a)) 

(3  2) 
X2{co)  =  S(coyj0}L  +N2(co), 

where  the  uppercase  letters  represent  Fourier  transforms  of  the  respective  signals  in  Eq. 
3.1. 

If  the  transient  arriving  at  sonabuoy  2  is  different  from  the  transient  arriving  at 
sonabuoy  1 ,  the  received  signals  are  given  by: 

xl(k)  =  s(k)+j]l(k) 

x2(k)  =  r(k-L)+rj2(k), 

where  r{k)  is  the  other  transient  arriving  at  sonabuoy  2  with  time  delay  L  relative  to  s(k). 
The  corresponding  expressions  in  the  frequency  domain  are 

xAco)=S(co)+N]{cd) 

(3  4) 
X2{co)=R(co)e-]a}i  +N2(co). 

If  there  is  no  transient  present  at  the  second  sonabuoy,  the  two  received  signals 


are 


xx{k)  =  s(k)+n,(k) 
x2(k)  =  7j2(k). 

In  this  case  the  signal  received  at  sonabuoy  2  consists  entirely  of  noise  and  the  frequency 
domain  representation  of  the  two  signals  is  thus: 

X2{o))=N2{co). 

B.         MOMENTS  AND  CUMULANTS  OF  A  RANDOM  PROCESS 

For  the  purposes  of  this  thesis,  the  following  notation  and  definitions  are  used. 
The  n"'  moment  of  a  real  stationary  random  process  is  [Ref  7:p.  15] 
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<  fa > Ti  >-•,*-, ,)  =  £{*(*)*(*  +  r,  )..jc{k  +  zn_x )},  (3.  7) 

where  £  denotes  statistical  expectation  and  ti  represents  the  i  lag.   Note  that  m\  (r)  is 
the  ordinary  autocorrelation  function  . 

In  the  analysis  of  signals  using  higher  order  statistics,  cumulants  rather  than 
moments  are  generally  used.  Cumulants  of  order  n  are  defined  by  certain  linear 
combinations  of  products  of  moments  of  order  n  and  lower  [Ref  7:p.  15].  Cumulants  of 
Gaussian  random  processes  have  the  distinction  that  they  are  all  zero  for  order  n  greater 
than  2.  Thus  signals  imbedded  in  additive  gaussian  noise  theoretically  appear  naked 
when  subjected  to  analysis  using  higher  order  cumulants.  The  n'h  order  cumulant  is 
denoted  by 

c" (r, , r2 , ,rB_, ,)  =  Cum[x{k)x(k  +  r,  )...jc(k  +  zn_x )]  (3.  8) 

For  a  more  detailed  definition  of  cumulants  the  reader  is  referred  to  [Ref  7:p.  9]. 

The  following  relationships  exist  between  moments  and  cumulants  for  stationary 
random  processes  [Ref  7:p.  9]: 

a.  First  Order. 

clx=m[=E{x(k)}  (3.9) 

b.  Second  Order. 

c2AT])  =  m;(T])-(m]x)2  (3.10) 

c.  Third  Order 

c] (r, , r2 )  =  m2x (r, ,r3)-(mlx)  [m] (r, )  +  m] (t2 )  +  m] (r2  - r, )]+  l{m\ ) 

(3.11) 
Using  these  relationships,  zero  mean,  white  Gaussian  noise  has  the  characteristics  listed 
in  Table  1 . 
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n 

n 

m 

n 
C 

1 

0 

0 

2 

a2 

a2 

3 

0 

0 

4 

3a4 

0 

Table  1. 


Moments  and  cumulants  for  white  gaussian  noise  for  r,  =  r,  =  r,  =  0, 


C. 


N-TH  ORDER  "MOMENTS"  OF  A  DETERMINISTIC  SIGNAL 


In  this  thesis,  we  consider  transients  to  be  deterministic  signals.  Consequently, 
concepts  such  as  statistical  moments  are  not  defined.  However  certain  operations  in  the 
time  domain,  which  are  analogous  to  estimating  moments  for  realizations  of  a  random 
process,  are  still  useful  for  deterministic  signals.  Some  authors  (e.g.,  Nikias  and 
Petropulu  [Ref  7])  have  referred  to  these  operations  as  computing  "moments  and  higher 
order  spectra  for  deterministic  signals."  Since  some  of  the  techniques  we  have  adapted 
are  due  to  authors  using  this  concept,  we  will  adopt  this  concept  here  as  well. 

In  general  the  n'h  order  moment  for  an  energy  signal,  x(k),  is  defined  as  [Ref  7:p. 
78] 


m 


;(r1,....v,)=   I  x(k)x(k  +  rXx{k  +  Tn-:\ 


*  =  -CC 


and  for  a  power  signal  [Ref  7:p.  100]  as 


I    J+N-\ 

<(rP--Vi)  =  -   £   x{k)x{k  +  zx)....x{k  +  zn_x), 

N    k=J 


(3.  12) 


(3.  13) 


where  J  is  an  arbitrary  starting  point  of  the  summation  and  ./Vis  the  period  of  the  signal. 
From  these  expressions,  moments  can  be  considered  as  a  measure  of  the  degree  of 
similarity  between  a  signal  and  delayed  or  advanced  replicas  of  itself.  The  n'h  order  cross 
moment  for  n  energy  signals  is  defined  as 


in 


"   ..., .(^-..r-i)-   5  *,  (*k (*  + r,  )...*„ (*  +  rn_,)  (3.  14) 


k=-°o 
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The  n'h  order  moment  spectrum  is  the  Fourier  transform  of  Eq.  3.12  and  can  be 
expressed  as  Fourier  transforms  of  the  signal.  For  energy  signals,  this  is  given  by  [Ref 
7:p.85] 

Mnx(o)1,...cDB_l)  =  X(o)iy.X(a>n_1)X\a)]  +  ..  +  »_»).  (3.  15) 

which  can  be  written  in  terms  of  magnitude  and  phase  as 

\Mnx  {co,  ,...*v,  )|  =|  X{a>x )  I ..  I  X{mn_x )  ||  X(co,  + ..  +  con_x )  | 

(3.  16) 
x¥"x(o),,...con_])  =  (p(col)  +  ...  +  <p(con_l)-^(a)]..,con_}), 

where    M  "(&>,,... &>,,_,)    is  the  magnitude  term  and   x¥"(o)i,...o)n_i)   is  the  phase  term. 

Orders  n=2,  3  and  4  are  important  special  cases  of  moments.  In  the  frequency  domain 
these  orders  have  been  termed  Power  Spectral  Density  (PSD),  Bispectrum,  and 
Trispectrum,  respectively. 

D.         SECOND-ORDER  MOMENTS 

1.  Definitions  of  Moments  and  Spectra 

Using  Eq  3.12,  the  autocorrelation  of  a  deterministic  signal,  x(k),  is  written  as 

ml(r)=  I  x(k)x(k  +  r).  (3.  17) 

k=— oo 

From  Eq  3.16,  the  magnitude  and  phase  components  of  the  PSD  are  given  by 

\M;U  =  \x(cof 

^2(a>)  =  0. 

From  these  expressions,  it  can  be  seen  that  the  PSD  is  an  even  function  with  no  phase 
information.  Similarly,  the  cross-correlation  of  two  deterministic  signals  is 


m;Jr)=  Z  X]{k)x2{k  +  r).  (3.19) 

The  corresponding  cross-spectrum  is  the  Fourier  transform  of  Eq  3.19: 

M*iX2(a>)  =  X1{cD)X2*(a>),  (3.20) 

which  can  be  written  in  terms  of  magnitude  and  phase 
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MlXco\=\Xx{c)}\X2(co\ 

'  (3.21) 

2.  Second  Order  Moments  of  the  Received  Signals 

Since  cross-correlation  is  used  extensively  throughout  the  thesis,  it  is  important  to 
determine  the  cross-correlations  for  the  three  signal  cases  presented  earlier.  Before  doing 
this,  let  us  first  investigate  the  cross-correlation  between  the  two  noise  sources  77,  (k)  and 
r\2  (k)  at  the  two  sonabuoys  of  interest.  The  cross-correlation  and  the  corresponding 
cross-spectrum  are  given  by: 

00 

M^  (a>)  =  %{mL  to}  =  N,  {co)N2  (co). 
The  expectation  of  this  term  is  zero  since  the  two  white  noise  sources  are  uncorrelated. 
However,  the  term,  M1    ,  itself  is  not  zero.    Applying  this  to  the  case  where  the  same 
transient  arrives  at  the  two  sonabuoys,  the  cross  spectrum  is 

MlxXco)  =  \S(cofe-J<oL  +M2mJ<*)+K  (*\       ■    (3-  23) 
where  M2    are  the  cross  terms.   From  this  equation,  it  can  be  seen  that  the  time  delay  L 

is  imbedded  in  the  linear  phase  term  e~JcoL .  This  linear  phase  term  is  therefore  important 
in  finding  the  TDOA  between  the  two  signals. 

For  the  case  where  different  transients  arrive  at  the  two  sonabuoys,  the  cross- 
spectrum  is 

Ml,  (*)  =  |sH|i?He«>^<«  V*"  +  Ml,  (co)  +  Ml  W         (3.  24) 

For  this  case  it  can  be  seen  that  there  is  the  same  linear  phase  term  e~ J0jL  as  that  of  Eq 
3.24.  Therefore  the  same  TDOA  result  would  be  observed  for  Eq  3.24  and  3.25  even 
though  the  transients  are  different.  The  linear  phase  can  therefore  not  be  used  as  a  means 
to  discriminate  between  two  transients.  In  the  last  case,  where  there  is  only  noise  present, 
the  cross-correlation  is 

M;Jco)  =  M;Jco).  (3.25) 
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For  this  case,  it  is  expected  that  the  phase  will  be  random,  and  therefore  no  linear  phase 
term  will  exist. 

E.         THIRD-ORDER  MOMENTS  -  BICORRELATION  AND  BISPECTRUM 

1.  Definition  of  Moments  and  Spectra 

The  third-order  moment  of  a  signal  is  called  the  bicorrelation.  From  Eq  3.12,  the 
auto-bicorrelation  of  a  deterministic  real  signal  is  defined  as 


ml (ri ' T2 ) =  £  x(k)x(k  +  r, )x(k  +  t2)  (3 .  26) 

and  from  Eq  3.14  the  cross  bicorrelation  is 

CO 

Kv,  (r> ' r2 )  =  I  *,  {k)x2  (k  +  r,  )x3  (k  +  t2  ).  (3.  27) 

k  =-oo 

For  this  thesis,  where  there  are  only  two  data  streams,  the  cross-bicorrelation  will  consist 
of  only  two  signals,  xx{k)  and  x2\k).  The  cross-bicorrelation  functions  of  interest  are 


<*,„  (*i >t2)  =   I  x,  (*)*,  (A:  +  r,  )x2 (*  +  r2 )  (3.  28) 

1    '    2  t=-oo 

and  the  terms  rnrr  ,  ,  ml  ,r   andm^ , ,  ,  which  are  defined  in  an  analogous  way. 
The  auto-bispectrum  is 

Mj  (©, ,  tf>2 )  =  *(a>,  )x(<y2  )jf  *  (fi>,  +  co2 ),  (3.  29) 

or  in  terms  of  magnitude  and  phase 

tyl  (®1  >*>2  |  =  |*(*>1  ||^(^2  )||^fe  +  ®2 ) 
y*(Ci>l><i>2)=<f>xM+4>XM-<f>x(<t>l+<»2)- 

As  can  be  seen,  by  comparing  Eq  3.18  and  Eq  3.30,  a  distinct  difference  between  the 


second  order  moment  spectrum  and  the  bispectrum  is  that  the  bispectrum  contains  phase 
information.  The  cross-bispectrum  is 

M\x,M  fo  '  ^2  )  =  X2  (P\  )X3  {®2  )X*fa+G>2\  (3.31) 

or  in  terms  of  magnitude  and  phase 
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M*m  fe  ,et>2U  \X2 (a,,  }\X3 {co2 }\X,  fa  +  co2 } 

(3.  32) 

^x,  fo  >  °>2  )  =  K  fo  )  +  #c,  fo  )  -  ^  fo  +  ^2  ) 

2.  Third  Order  Moments  of  the  Received  Signals 

Both  second  order  moments  and  third  order  moments  are  used  extensively 
throughout  this  thesis.  It  is  therefore  important  to  develop  the  third  order  moments  for 
the  three  signal  cases.  To  provide  a  motivation  for  the  algorithms  to  be  used,  let  us 
consider  a  situation  in  which  the  noise  is  identically  zero.  In  reality  the  noise  terms 
TV, (<i>)  and  N2(a>)  are  not  zero,  although  the  expected  value  of  their  higher  order 
moments  defined  earlier  vanish  when  the  noise  is  Gaussian. 

For  the  case  where  the  same  transient  arrives  at  both  sonabuoys,  the  signals  in  the 
frequency  domain  are  given  by  Eq.  3.2.  Using  Eq.  3.32,  the  bispectrum  of  these  signals 
is  given  by 

Wm  K  >  ®2  }  =  \S(°>i  l^-  ]\S(o)x  +  co2 } 

(3.  33) 

^ivi  fo ,co2)  =  0s{co])-coL  +  <ps (co2 ) - <ps (fi>,  +  co2 ), 
where  (f)s  is  the  phase  of  the  signal  S(a>).   Equation  3.33  is  the  expression  for  the  cross- 
bispectrum  between  the  data  streams  at  sonabuoy  1  and  sonabuoy  2  in  the  absence  of 
noise.     The  signal  phase  terms,   (j)s (a)] )  +  (f)s (co2 ) - (f)s (co]  +co2),  can  be  eliminated  by 
making  use  of  the  auto-bispectrum  at  sonabuoy  1,  which  is  given  by 

|M;  fo  ,cd2  )|  =  \Sfa  )\S{co2  pico,  +a>2] 

^  fo  *<*>l)  =  4sM+  <f>S  (*>2  )-  ^S  (®1  +  ^2  ) 

This  expression  is  then  used  to  normalize  the  cross-bispectrum  and  results  in  the  identity: 

zfa,fl>2)=  ,;;;      :  =e-^L.  (3.35) 

From  Eq  3.35  it  can  be  seen  that  by  normalizing  the  bispectrum  all  the  phase  terms  are 
cancelled  except  for  the  linear  phase  term,  e'Jto'L  . 
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For  the  case  of  different  transient  arrivals,  the  cross-bispectrum  can  be  obtained 
from  Eq.  3.4  and  Eq.  3.32  as  follows 

M"Wl  (*>. > *>2  ]  =  \R(^  %S(a>2 pico,  +o)2) 

^XM  (°>l  ><°2  )  =  fa  (®1 )  -  °>L  +  0S  (&2  )  -  <t>S  («1   +  <°2  \ 

where  fa{&)  is  the  phase  of  the  signal  R(co).  In  this  case  the  normalized  cross- 
bispectrum  takes  the  form 

Note  that,  unlike  the  first  case  (see  Eq.  3.35),  the  linear  phase  term  cannot  be  separated 
from  the  other  phase  terms.  In  both  cases  however,  the  two-dimensional  bispectrum  is 
reduced  to  a  function  of  a  single  variable  a>l . 

F.         SIGNAL  PROCESSING  ALGORITHMS 

1.  Bispectrum  Linear  Phase  Detector 

The  analysis  of  the  previous  section  shows  that  when  the  same  signal  is  arriving 
at  both  sonabuoys,  a  linear  phase  term  is  present  in  the  normalized  function  l(co] ,co2). 
The  time  delay  L  can  now  be  extracted  using  the  following  ad  hoc  method  developed  for 
time  delay  estimation  [Ref  16]. 

To  estimate  the  delay  L,  a  third-order  "hologram"  transformation  is  required. 
This  is  defined  to  be  [Ref  7:p.  324] 


JL       Ji. 

T(t)=  J  jl(oliQ)2)e~j<0,TdQ)ld(D2.  (3.  38) 


Since  /  reduces  the  bispectrum  to  one  dimension,  the  third-order  hologram  is  a  one- 
dimensional  Fourier  transform  over  the  dimension  containing  the  linear  phase  term, 
followed  by  an  integration  over  the  second  dimension.  Note  that  when  the  same  signal  is 
present  at  both  sonabuoys,  Eq.  3.38  takes  the  form 
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T(t)=  J  \ejloAL-T)dco,dco2.  (3.39) 


Since  all  the  elements  of  this  second  dimension  are  in  phase  they  add  up  constructively, 
giving  a  strong  peak  at  T  =  L.  Since  we  are  working  in  the  discrete-time  domain  the 
third-order  hologram  can  be  rewritten  as 

M-\  M-\ 


^W=SS/(®,^2>"Mr.  (3-40) 


(O,  =0<i>T  =0 


The  absolute  value  of  T(z)  will  display  a  strong  peak  at  the  location  of  the  of  the  time 
delay  between  the  two  signals  of  sonabuoy  1  and  sonabuoy  2. 

For  the  case  where  the  two  transient  arrivals  are  different,   the  third-order 


hologram  takes  the  form: 


T{T)=  Y^^ei^^)e^M^^d<old<D1.  (3.  41) 


In  this  case,  the  third  order  hologram  contains  extra  phase  terms  fa  and  (f>s,  which  will 
either  add  in  phase  or  out  of  phase.  Thus,  in  general,  T(t)  will  not  exhibit  a  strong  peak. 

2.  Bispectrum  Discriminator 

The  bispectrum  linear  phase  detector  can  also  be  used  as  a  discriminator  by 
applying  a  simple  threshold  technique  to  the  third  order  hologram.  The  threshold 
technique  is  applied  by  first  noticing,  from  Eq  3.41,  that  if  the  transients  are  different  the 

\R(o){ } 
hologram  contains  phase  terms  (f)R-  fc  -coL  and  magnitudes  \— -. — ^ .  On  the  other  hand,  if 

|%J 
the  transient  signals  are  the  same,  the  hologram  contains  only  the  linear  phase  term  caL 
(Eq  3.39).  The  phase  and  magnitude  terms  of  Eq  3.41  will  add  constructively  or 
destructively  to  produce  peaks  and  valleys  to  the  bispectral  hologram.  If  the  signals  are 
the  same,  there  will  be  only  one  peak  in  the  hologram,  and  that  peak  will  be  at  the  delay 
L. 

In  summary,  one  can  expect  that  if  the  SNR  is  sufficiently  large  and  the  signals 
are  the  same,  a  peak  at  delay  L  will  dominate  the  hologram.  On  the  other  hand  if  the 
signals  are  different,  there  will  be  a  maximum  value  at  delay  L  as  well  as  other  extrema  at 


delays  other  than  L.    These  other  extrema  will  be  larger  than  the  noise  and  can  therefore 
be  used  to  discriminate.  A  discrimination  algorithm  is  therefore  proposed  as  follows: 


Step  1: 

Estimate  T(t) 

Step  2: 

Find  the  maximum  value  of  the  hologram 

£>  =  argmax  T(t) 

t      \ 

(Q  is  the  magnitude  of  the  hologram  at  the  estimated  TDOA). 

Step  3: 

If  r0  is  the  value  of  r  that  produces  the  maximum  in  step  2  then  compute 

Q2  =  argmax  T(t)-Q 

T  V 

Step  4: 

Define  a  threshold  using  Eqs.  3.42  and  3.43  as 


1 


(n-\ 


Qt=^—yj^)-q-q,« 


\i=o 


N-2 
Step  5: 
Compare  QT  to  the  difference  Q  -Q2.  If 

Qt>Q-Qi 
the  signals  are  considered  to  be  different;  if 

QT<Q-Q2 

the  two  signals  are  considered  to  be  the  same. 


(3.  42) 


(3.  43) 


(3.  44) 


(3.  45) 
(3.  46) 
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IV.    SUBSPACE-BASED  SIGNAL  PROCESSING 

Signal  subspace  methods  have  been  used  in  a  variety  of  problems  for  estimating 
spectral  lines  and  direction  of  arrival  in  sensor  arrays.  A  good  introduction  to  these 
methods  for  spectrum  estimation  can  be  found  in  [Ref  17].  Recently,  the  use  of  subspace 
methods  has  been  explored  for  estimating  TDOA  [Ref  20],  [Ref  21].  This  thesis  uses  the 
same  approach  as  [Ref  20]  for  TDOA  estimation  and  further  adapts  the  subspace  method 
for  transient  discrimination. 

The  discussion  begins  by  defining  the  subspace  model  for  the  TDOA  problem  and 
then  moves  on  to  the  MUSIC  method,  which  was  used  for  delay  estimation  and 
discrimination.  The  chapter  finishes  with  a  discussion  of  how  the  MUSIC  method  was 
adapted  and  used  for  discrimination 

A.         SUBSPACE 

To  understand  the  principle  of  subspace  methods  and  their  application  to  the  topic 
of  this  thesis,  consider  the  case  where  the  same  signal  is  arriving  at  two  separate 
sonabuoys  (see  Eq.  3.1).  The  frequency  domain  representation  of  Eq.  3.1  is  given  by  Eq 
3.2.  From  Eq.  3.2,  it  can  be  seen  that  the  linear  phase  term  eja)L  contains  the  delay  L 
between  the  transient  arrivals  at  sonabuoy  1  and  sonabuoy  2.  The  cross-spectrum  for  the 
two  received  signals,  given  by  Eq  3.23,  is  repeated  here  for  convenience 

M2X]Ja>)  =  \S(cofe-^  +M;Jco)  +  M2xr}  {co).  (4.  1) 

It  can  be  argued  [see  Ref  20]  that  when  Eq  4.1  is  sampled  at  equally-spaced  values  in 
frequency,    the   resulting    data    sequence,    y(k)  =  M\x  (cok ) ,    satisfies    the    conditions 

necessary  to  apply  a  signal  subspace  model  [Ref  17].  Specifically,  an  NxN  correlation 
matrix  is  estimated  for  the  data  y(k),  and  the  corresponding  iV-dimensional  vector  space  is 
divided  into  signal  and  noise  subspaces.  The  approach  followed  to  estimate  TDOA  using 
subspace  methods  is  to  project  a  steering  vector  of  the  form 
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d(/)= 


(4.2) 


,jl(N-\) 

onto  the  noise  subspace  and  plot  the  result  as  the  linear  phase,  /,  is  varied.  When  the 
parameter  /  is  equal  to  the  true  delay  L  between  the  sonabuoys,  the  projection  onto  the 
noise  space  is  zero.  In  subspace  algorithms,  such  as  MUSIC,  this  null  projection  can  be 
used  to  estimate  TDOA. 

Unfortunately,  the  cross-spectrum  for  the  case  where  different  transients  are 
arriving  at  sonabuoys  1  and  2  also  contains  a  linear  phase  term  (see  Eq.  3.24).  Therefore, 
the  basic  subspace  techniques  are  not  useful  for  transient  discrimination  without  suitable 
modification. 


B. 


MUSIC 


In  the  application  of  subspace  techniques  to  linear  phase  detection  between  two 
transients,  the  MUSIC  (Multiple  Signal  Classification)  method  developed  by  Schmidt 
[Ref  18]  is  used.  A  brief  outline  of  the  method  is  provided  here. 

The  MUSIC  algorithm  is  implemented  in  the  frequency  domain  by  first  obtaining 
the  cross-spectrum  of  the  two  signals  received  at  the  sonabuoys.  Since  M\x  {co)  is 
formed  using  a  DFT,  this  data  is  available  at  samples  cok ,  k=  0,1,2. ..Ndft-1,  where  Ndft 
is  the  size  of  the  DFT.  From  this  frequency  domain  data,  a  data  matrix  M^  is  formed, 
and  the  correlation  matrix  is  estimated  as 

This  can  be  decomposed  into  an  orthonormal  eigenvector  matrix 

E  =  [E,gE_J  (4.4) 

and  a  diagonal  eigenvalue  matrix 


See  [Ref  17]  for  various  methods  of  forming  a  data  matrix 
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A  = 


A^0 


(4.5) 


.0      Ano„ 

The  columns  of  Es,o  and  EB0«e  fonn  an  orthogonal  basis  set  and  define  the  signal  and 
noise  subspaces,  respectively.  It  is  therefore  possible  to  write  Rx  as 

R    -  EAE*r  =  E  .  A  .  E*T  +  E    -A    .  E*r    .  (4.  6) 

X         •-"■*-"  sig       sig       sig  noise       noise      noise  \  f 

The  columns  of  the  eigenvector  matrix  can  be  used  to  form  projection  matrices  for  the 
signal  and  noise  subspaces  of  the  form  [Ref  17:p.  623]: 


s<g 


E  .  E*r 

sig      sig 


P       =E      E*T.  . 

noise  noise      noise 


(4.7) 


When  the  projection  matrix  P  is  multiplied  by  a  vector  d,  the  result  d  =  Pd  is  the 
projection  of  d  into  the  corresponding  subspace  (see  Figure  5). 


Signal 
Basis  Vector 


Noise 

Basis  Vector  2 


Noise 

Basis  Vector  1 


Figure  5.         Projection  of  vector  d  onto  noise  subspace  by  MUSIC. 


As  discussed  previously,  subspace  methods  make  use  of  the  fact  that  the 
projection  of  the  steering  vector  d(/)  ,  given  by  Eq.  4.2  onto  the  noise  subspace  is  zero. 
The  MUSIC  algorithm  [Ref  17:p.  627],  in  particular,  evaluates  the  quantity 

1  1  1 


MU 


d*r(/)Pno„ed(/)      d*r(/^ee;rd(/)' 


(4.8) 
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where  PMU  is  the  MUSIC  indicator  function  for  d(/)and  e2...eN  are  the  eigenvetors 
spanning  the  noise  subspace  (the  eigenvector  e,  spans  the  signal-subspace  which  is  one- 
dimensional)  The  value  of  /  where  PMU  exhibits  a  sharp  peak  determines  the  TDOA 
between  the  signals. 

C.         SUBSPACE  DISCRIMINATOR 

A  subspace  discriminator  was  developed  by  comparing  projections  onto  multiple 
subspaces  using  the  MUSIC  algorithm.  The  following  considerations  are  proposed  to 
motivate  this  discriminator. 

Since  the  subspace  method  in  this  thesis  is  based  on  the  cross-spectrum  between 
two  sonabuoys  (see  Eq.  4.1)  it  is  important  to  investigate  the  differences  between  the 
cross-spectrum  for  the  case  where  the  same  transient  arrives  at  two  separate  buoys  and 
that  for  the  case  where  different  transients  arrive  at  two  separate  buoys.  This  will  be 
helpful  in  understanding  the  subspace  discriminator.  From  Eq.  3.24,  it  can  be  seen  that  if 
the  transients  are  different,  the  cross-spectrum  contains  phase  terms  (f>s{co)-tpR{co)-coL 

and  magnitudes  \R{cox  J\s(col  )j .  However,  if  the  transient  signals  are  the  same,  the  cross- 
spectrum  contains  only  the  linear  phase  term  coL  (Eq.  3.23).  Now,  if  the  subspaces 
formed  using  Eq.  3.23  are  compared  to  those  formed  using  Eq.  3.24  the  two  subspaces 
would  typically  be  rotated  with  respect  to  each  other  as  illustrated  in  Figure  6.  The 
implication  of  the  rotation  of  the  subspaces  due  to  the  phase  terms  of  Eq.  3.24  is  that 
when  the  parameter  /  of  the  steering  vector,  d(/),  is  equal  to  the  delay  L  between  the 
sonabuoys,  the  projection  onto  the  noise  subspace  will  not  be  zero  as  in  the  case  when  the 
same  transient  arrives  at  the  two  sonabuoys.  By  comparing  projections,  it  would  be 
found  that  the  differences  in  magnitude  of  the  projections  would  be  large  if  the  transients 
were  the  same.  On  the  other  hand  if  the  transients  were  different,  the  differences  in 
magnitude  between  the  projections  would  be  small.  These  differences  in  magnitude  are 
used  to  discrimnate  between  transients. 

This  method  is  developed  further  by  noticing  that  the  noise  projection  matrix, 
P noise,  can  be  written  as: 
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Pn^-Ze<er=I-e,e;r,  (4.9) 

1=2 

where  I  is  an  NxN  identity  matrix. Let  us  further  define  Pa  as  the  projection  operator  for 
the  subspace  formed  by  leaving  out  the  k    basis  vector.  Thus 

P,  =2>,e*r=I-e,e7.  (4.10) 

M 

i*k 

Further  we  define  the  following  set  of  additional  indicator  functions 

d    (l)Vkd(l)     d    (Hl-ekek  ]d(/j 
The  proposed  subspace-based  discrimination  algorithm  is  as  follows: 


Step  1 

Estimate  the  correlation  matrix  Rx  and  find  the  eigenvectors  e,. 

Step  2 

Compute  the  maximum  estimate 

( *\ 
Q  =  arg  max  P  .  (4.  12) 

*,/     V   ) 
Step  3 
If  ko  is  the  value  of  k  that  produces  the  maximum  in  step  2  then  compute 

(4.  13) 


Q2 

=  arg  max 

k,i      V 

k*k0 

p  . 
J 

Step  4 

Compare  the  difference  Q 

-Qz- 

If 

Q- 

Q2 

<i 

the  signals 

are  considered  to  be  different. 

if 

Q- 

Q2 

<i 

the  signals  are  considered  to  be  the  same. 

(4.  14) 
(4.  15) 
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A    Signal 
Basis  Vector 
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Basis  Vector' 


Noise 

Basis  Vector  1" 


Noise 

Basis  Vector  1 


Noise 

Basis  Vector  2 


Noise 

Basis  Vector  2' 


Figure  6.  Relative  rotation  of  subspace  due  to  phase  terms  of  Eq.  3.24. 
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V.    SIMULATION  RESULTS 


A.         SIMULATION  CONDITIONS 


1.  Synthetic  Transients 


In  this  thesis  a  set  of  four  synthetic  transients  was  used  to  evaluate  each  of  the 
algorithms  described  in  the  previous  chapters.  The  four  transients  were  a  CW  sinusoidal 
pulse,  an  exponentially  decaying  sinusoidal  pulse,  a  linear  phase  modulated  (LFM)  pulse 
and  a  simulated  finback  whale  transient. 

Each  transient  can  be  described  as  an  amplitude-  and  phase-modulated  sinusoid, 
given  by  the  general  expression  [Ref  19:p.  202] 

p(t)  =  a(t)  cos  0(t),  (5.1) 

where  a(t)  is  the  time- varying  amplitude  or  "envelope"  of  the  signal  and  0(t)  is  the  time- 
varying  angle.  The  angle  Oft)  is  of  the  form  0(t)  -  27tfct  +  y{t)  so  that 

p{t)  =  a{t)cos{(oct+y(t)),  (5.2) 

where  fc  is  the  center  frequency  and  yft)  is  the  phase  modulation.  Each  of  these  transients 

is  described  in  more  detail  below. 

For  the  CW  pulse,  the  envelope  and  phase  modulation  are  given  by: 

fl,    0  <  t  <  T 
a(t)  =  \  (5.3) 

[0,     otherwise 

y(t)  =  0  (5.4) 

where  T  is  the  pulse  length.    A  plot  of  this  transient  is  given  in  Figure  7  for  T  =  40  ms 

and  fc  =50 Hz. 
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Figure  7.  CW  Pulse. 

The  exponentially  decaying  sinusoidal  transient  is  described  by  the  following 
envelope  and  phase  modulation  terms 

,x    IV",   o<r<r 

a(t)=      A  L  •  (5-5) 

[0,  otherwise 

K0  =  0,  (5.6) 

where  T  is  the  pulse  duration  and  a  is  the  attenuation  constant.  This  transient  is  shown  in 
Figure  8  for  a  =  200  and  T  =  40  ms . 
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Figure  8.  Exponentially  Decaying  Sinusoidal  Transient. 


For  the  Linear  Phase  Modulated  (LFM)  transient,  a  Gaussian  envelope  of  the 


form 
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'(')= 


e    /2°2 ,  0  <  t  <  T 
0,  otherwise 


(5.7) 


was  used.  The  phase  is  a  quadratic  function  of  time: 

In-Af-r 


7(0  = 


27 


(5.8) 


where  Af  is  the  desired  change  in  instantaneous  frequency  over  the  interval  T.  A  plot  of 


.2        T 


this  transient  for  T  -  40 ms  , cr"  =  —and  A/  =  100#z  is  shown  in  Figure  9. 
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Figure  9.  LFM  Pulse. 

The  finback  whale  transient  is  the  most  complicated  of  all  the  transients  that  were 
synthesized.  This  transient  is  modeled  as  follows  [Ref  8]: 

\3/T)-t,         0<t<T/3 
1,  773<r<2773 

3-(3/T)-t,     2T/3<t<T 
0,  otherwise 

y(t)  =  2x-23.eh{]*:,23:)'  (5.10) 

with  fc  -  0.  This  transient  is  plotted  in  Figure  10  for  T  -\  s  .  The  instantaneous 
frequency  decreases  nonlinearly  from  23  Hz  to  18  Hz  with  the  most  rapid  decrease 
occurring  initially  and  the  rate  of  decrease  becoming  smaller  with  time. 
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Figure  10.        Finback  Whale  Transient. 


2. 


Sisnal-to-Noise  Ratio 


During  the  evaluation  of  the  algorithms,  the  synthetic  signals  described  above 
were  subject  to  additive  white  Gaussian  noise.  The  signal-to-noise  ratio  (SNR)  is  defined 
as  follows.  Let  Ptmnsienl  represent  the  average  power  in  the  signal: 


N  -1 


1       'v,-i 

P  =  —  Y\x(n)\2 

1  transient  ,  r       Lu\     \     >  \     ' 


(5.11) 


Note  that  the  transient  signal  power  is  normalized  by  the  length  of  the  transient  Ns  and 
not  the  entire  observation  time.  This  is  done  to  prevent  the  SNR  from  changing  with 
changing  observation  time.  The  SNR  in  dB  is  then  defined  as 


SNR  =  101ogI0  P'rans,en'  =  lOlog.o  P'~l  dB 


P 


(5.  12) 


Noise 


where  <J2noise  is  the  noise  variance  used  in  the  generation  of  the  Gaussian  white  noise. 


B. 


RESULTS 


To  evaluate  the  algorithms,  two  sets  of  experiments  were  used,  one  to  perform 
TDOA  estimation  and  the  other  to  perform  discrimination.  For  the  first  set  of 
experiments,  the  same  transient  arrives  at  both  sonabuoys.  For  this  case  all  the  transients 
were  evaluated  in  turn,  over  a  range  of  SNR  values.   For  the  second  set  of  experiments, 
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all  combinations  of  transient  arrivals  at  sonabuoys  1  and  sonabuoys  2  were  tested  for 
discrimination. 

In  all  of  the  experiments,  the  transients  were  set  to  the  following  lengths:  20  ms 
for  the  CW  sinusoidal  pulse,  40  ms  for  the  exponential  decaying  sinusoid  and  LFM  pulse, 
and  1  s  for  the  whale  transient.  The  delay,  L,  between  the  transients  was  always  kept  at 
40  samples,  which  is  equivalent  to  200ms.  A  sampling  rate  of  four  times  the  required 
Nyquist  frequency  was  used  for  each  transient.  This  satisfies  the  requirement  for  the 
bispectrum,  which  requires  a  sampling  rate  of  three  times  the  maximum  frequency  [Ref 
8].  The  total  length  of  each  signal  used  was  256  samples  or  1.28  s.  The  parameters  are 
summarized  in  Table  2. 


/c(Hz) 

T 

L  (delay) 

A/ 

a 

CW  Pulse 

50 

20ms 

200ms 

Exponential 
Decaying  Sinusoid 

50 

40ms 

200ms 

200 

LFM 

50 

40ms 

200ms 

1kHz 

Whale 

Is 

200ms 

Table  2. 


Transient  parameters  used  in  experiments. 


1. 


TDOA  Estimation 


The  TDOA  results  will  be  discussed  as  follows.  First,  an  example  of  the  basic 
results  using  the  exponentially  decaying  sinusoid  transient  will  be  given  for  both  the 
subspace  method  and  the  bispectrum  linear  phase  detector.  After  this,  certain 
"probability"  curves  for  TDOA  will  be  defined  and  presented  for  each  transient  starting 
with  the  subspace  method  and  ending  with  the  bispectrum  results. 

The  subspace  method  uses  the  MUSIC  algorithm  with  a  correlation  matrix  (Eq. 

4.3)  of  size  TV  =  60.    Figure  1 1  shows  a  plot  of  the  function  Pmu   (see  Eq.  4.8)  as  a 
function  of  /  using  a  SNR  of  12  dB.  Figure  12  is  a  plot  of  the  results  for  a  SNR  of  5  dB. 
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Figure  11.       Subspace  TDOA  estimation  for  an  exponentially  decaying  sinusoid 
using  an  SNR  of  12  dB  and  correlation  matrix  size  of  N=60. 

As  can  be  seen  from  these  figures,  the  algorithm  estimates  the  correct  TDOA  of  40 
samples  or  200  ms  for  a  SNR  of  12  dB,  and  yields  a  completely  incorrect  value  of -63 
samples  for  an  SNR  of  5  dB. 
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Figure  12.        Subspace  TDOA  estimation  for  an  exponentially  decaying  sinusoid 
using  an  SNR  of  5  dB  and  correlation  matrix  size  of  N=60. 

The  rationale  for  using  a  large  correlation  matrix  is  twofold.  First,  if  the 
correlation  matrix  is  too  small,  the  TDOA  estimate  tends  to  be  less  accurate.  A  typical 
result  is  shown  in  Figure  13  for  the  12  dB-SNR  case  using  a  correlation  matrix  of  size  N 
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=  5.  As  can  be  seen  in  this  figure,  a  TDOA  of  38  samples  (instead  of  the  correct  value  40 
samples)  was  estimated,  and  the  peak  is  broader  than  that  shown  in  Figure  1 1 .  Secondly, 
for  the  subspace  method  to  function  as  a  discriminator,  it  was  found  that  better  results 
were  achieved  at  lower  SNR,  when  the  correlation  matrix  and  thus  the  observation  space 
was  larger. 
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Figure  13.       Subspace  TDOA  estimation  for  an  exponentially  decaying  sinusoid 
using  an  SNR  of  12  dB  and  correlation  matrix  size  of  N  =  5. 

For  the  bispectrum  linear  phase  detector,  no  windowing  was  applied  to  the  data. 
Windowing  is  usually  applied  to  obtain  smooth  estimates  of  the  bispectrum  [Ref  7:p. 
126].  Due  to  the  short  data  lengths  used  in  this  thesis,  however,  it  was  found  that 
windowing  did  not  improve  the  results  and  was  therefore  not  used.  Typical  results  for  the 
bispectrum  linear  phase  detector  are  shown  in  Figure  14  and  Figure  15  for  the 
exponentially  decaying  sinusoid  using  SNR  values  of  12  dB  and  5  dB.  In  both  cases  the 
TDOA  is  indicated  by  the  maximum  value  of  the  hologram,  which  in  turn  is  a  single 
sample.  The  performance  is  similar  to  that  of  the  subspace  method.  At  12  dB  the  result 
is  exactly  correct  while  at  5  dB  the  method  gives  completely  erroneous  results. 
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Figure  14.        Bispectrum  TDOA  estimation  for  an  exponentially  decaying  transient 

using  an  SNR  of  12  dB. 
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Figure  15.        Bispectrum  TDOA  estimation  for  an  exponentially  decaying  transient 

using  an  SNR  of  5  dB. 

Figure  16  shows  typical  TDOA  results  using  the  classical  generalized  cross- 
correlation  methods  (GCC)  [Ref  2]  at  a  SNR  value  of  12dB.  ROTH  [Ref  13],  SCOT  [Ref 
14]  and  PHAT  [Ref  15]  weighting  was  used  in  addition  to  straight  cross-correlation  (no 
special  weighting  was  used).  As  can  be  seen  in  Figure  16  the  correlation,  SCOT  and 
PHAT  algorithms  give  strong  peaks  at  the  estimated  TDOA  while  ROTH  gives  a  very 
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noisy  signal  with  no  dominant  strong  peak.  In  general,  it  was  found  that  these  techniques 
gave  less  accurate  TDOA  estimates  than  the  bispectral  and  subspace  algorithms,  even  at 
high  SNRs. 


50 

Estimated 

40 

TDOA  =  89 

g>30 

20 

'it  I 

,  iLullt 

i 

i  !i  1  i 

10 

mmM 

J 

w 

-100     -50 


0 

lag 


50       100 


0.25 

0.2 
0.15 

0.1 
0.05  \ 


Estimated 
TDOA  =  89 


ii 


III 


'MW 


!  ! 


i   Ii  1 1 


ni 


-100     -50         0        50       100 
lag 


0.4 
0.3 
0.2 
0.1  J 


Estimated 
TDOA;  4  82 


0.25 

0.2 
0.15 

0.1 
0.05  I 


Estimated 
TDOA  =  89 


u 


|M| 


-100     -50         0        50       100 
lag 


L 


II 


'*.  lIMi 


-100     -50         0        50       100 
lag 


Figure  16.        TDOA  estimates  using  (a)  cross  correlation,  (b)  ROTH,  (c)  SCOT  and 

(d)  PHAT  with  an  SNR  of  12  dB. 

TDOA-estimate  performance  of  the  subspace  and  bispectral  algorithms  was 
evaluated  using  simulations  of  one  thousand  realizations  at  a  fixed  SNR.  For  each 
realization,  a  different  white  noise  source  was  added  to  the  original  signal.  The  output  of 
these  simulations  was  used  to  create  a  statistical  measure  of  performance  for  TDOA, 
which  we  call  the  probability  of  correct  TDOA  Pj.  Pt  represents  the  fraction  of  1 000 
trials  for  which  the  algorithm  estimates  the  correct  time  delay  between  transients  received 
at  the  two  sonabuoys.  Pt  curves  were  generated  for  SNR  values  of  8  dB  to  20  dB  for 
both  the  subspace  method  and  the  bispectrum  linear  phase  detector. 

Figure  1 7  and  Figure  1 8  (a)  shows  the  Pj  results  of  the  subspace  algorithm  for  the 
CW  pulse.  As  can  be  seen  from  the  figures,  the  subspace  algorithm  gives  excellent 
results  for  SNR  >  17  dB,  where  a  PT  of  1 .0  is  achieved. 
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Figure  17.        Pt  versus  SNR  for  CW  pulse  using  the  subspace  linear  phase  detector. 

The  Pt  curves  for  the  exponentially  decaying  transient  and  the  LFM  transient  are 
shown  in  Figure  18  (b)  and  (c).  Once  again,  a  Pt  value  of  1.0  is  achieved  for  SNR  > 
17dB.  The  Pt  for  the  whale  transient  (Figure  18  (d))  shows  poorer  results,  achieving  a 
maximum  Pt  of  1.0  at  a  SNR  of  20  dB.  This  poorer  performance  of  the  algorithm  may 
be  be  attributed  to  two  factors.  First,  the  whale  transient  was  the  longest  of  all  the 
transients  used  (being  1  s  long)  thus  filling  most  of  the  data  observation  window. 
Secondly,  it  is  a  non-linear  function  giving  it  a  complicated  phase  structure.  It  is  also 
important  to  note  that  the  results  in  Figure  1 8  (d)  are  based  on  a  longer  data  observation 
length  than  used  for  the  other  TDOA  results  (see  Table  3). 

Figure  18  (a)  and  Figure  19  show  the  Pt  for  the  CW  pulse  using  the  bispectrum 
linear  phase  detector.  From  these  figures,  it  can  be  seen  that  a  PTof  1.0  is  never  achieved 
but  that  a  Pt  of  0.9  is  achieved  at  a  SNR  of  15  dB.  The  remaining  plots  in  Figure  18  (b), 
(c)  and  (d)  also  show  that  a  Pt  of  1.0  is  not  achieved  for  the  bispectrum  method.  For  the 
exponentially  decaying  sinusoid  (Figure  18  (b))  a  Pt  of  0.9  is  achieved  at  1 1  dB  while  for 
the  LFM  pulse  (see  Figure  18  (c))  a  Pt  of  0.9  is  achieved  at  12  dB.  As  in  the  case  of  the 
subspace  methods,  the  Pt  of  the  whale  transient  was  obtained  using  a  longer  data 
observation  window,  where  a  maximum  Pt  of  0.9  was  achieved  at  an  SNR  of  16  dB 
(Figure  18(d)). 
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Figure  18.        Combined  Pt  for  subspace  and  bispectrum  linear  phase  detectors:  (a) 
CW  pulse,  (b)  exponential  decaying  sinusoidal  transient,  (c)  LFM  pulse  and  (d) 

whale  transient. 
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Figure  19.        Pt  for  the  CW  pulse  using  the  bispectrum  linear  phase  detector. 
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Figure  1 8  gives  a  summary  of  the  TDOA  results  for  all  of  the  experiments.  As 
can  be  seen  in  Figure  18,  the  bispectrum  linear  phase  detector  has  higher  PT  values  at  low 
SNR  than  the  subspace  linear  phase  detector  for  the  exponential,  LFM  and  whale 
transients.  However  at  higher  values  of  SNR,  the  subspace  algorithm  gives  better  Pj 
results  (i.e.,  Pt  of  1.0)  which  is  never  achieved  using  the  bispectrum  linear  phase 
detector. 

In  comparing  the  two  methods  (subspace  and  bispectrum)  one  may  note  that 
although  the  bispectrum  method  never  achieves  a  Pt  of  1 .0,  in  three  out  of  the  four  cases 
it  reaches  the  value  of  Pt  =  0.9  earlier  than  the  subspace  method.  Table  3  compares  these 
results  by  listing  the  minimum  SNR  values  for  which  a  Pj  of  0.9  is  achieved  in  each  of 
the  test  cases. 


CW  Pulse 

Exponential 

LFM  Pulse 

Whale 

Subspace 

13dB 

15dB 

14dB 

18dB 

Bispectrum 

15dB 

lldB 

12dB 

16dB 

Table  3.  Minimum  SNR  required  to  achieve  PT  =  0.9. 

2.  Discrimination  Experiments 

The  transient  discrimination  results  will  be  discussed  in  a  fashion  similar  to  that 
of  the  TDOA  results.  The  discussion  will  therefore  start  by  giving  a  brief  introduction  of 
how  the  algorithms  were  implemented  by  using  the  CW  and  LFM  pulses  as  examples. 
After  this,  the  probability  curves  for  discrimination  will  be  discussed  for  every 
combination  of  transients. 

As  discussed  previously,  the  subspace  discriminator  makes  use  of  projections  into 

multiple  subspaces  and  compares  the  maximum  values  of  the  indicator  functions  Pk   Eq. 

A 

4.11  with  the  maximum  value  of  the  function  PMU    (Eq.  4.8).    It  is  found  that  if  the 

transients  arriving  at  the  sonabuoys  are  the  same  the  difference  in  peak  values  will  be 
large  as  shown  in  Figure  20(a).  If  the  difference  is  found  to  be  small  (i.e.,  <  1)  then  the 
two  transients  are  different,  as  shown  in  Figure  20(b).  Thus  in  making  these 
comparisons,  it  is  possible  to  discriminate  between  transients. 
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Figure  20.       Typical  Pk   estimates,  (a)  when  the  signals  are  the  same  and  (b)  when 
the  signals  are  different. 

In  the  case  of  the  bispectrum  discriminator,  the  phase  terms  fa  -  fa  (see  Eq  3.37), 
which  are  only  present  in  the  case  of  different  transient  arrivals,  add  up  constructively  or 
destructively  thus  giving  features  to  the  hologram  of  the  bispectral  linear  phase  detector. 
For  the  case  where  the  transients  are  the  same,  the  phase  terms  fa  -  fa  are  zero;  therefore 
no  extra  features  are  added  to  the  hologram.  If  the  difference  in  magnitude  between  these 
features  is  large  compared  to  the  average  of  the  hologram,  it  is  possible  to  discriminate 
between  transients.  These  characteristics  can  be  clearly  seen  in  Figure  21.  Figure  21(a) 
shows  the  single  peak  that  appears  for  similar  transient  arrivals  while  Figure  21(b)  shows 
the  multiple  peaks  and  extra  features  that  are  present  due  to  the  phase  terms  fa  and  fa  for 
different  transient  arrivals. 
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(a) 


Figure  21.       Bispectrum  linear  phase  detector.  Stem  plot  of  \t(tJ  for  (a)  similar 
transient  arrivals  and  (b)  different  transient  arrivals. 

The  experiments  for  transient  discrimination  were  as  follows.  Each  algorithm 
was  evaluated  using  simulations  of  one  thousand  realizations  at  a  fixed  SNR.  For  each 
realization,  a  different  white  noise  source  was  added  to  the  original  signal.  The  output  of 
these  simulations  was  used  to  estimate  the  probability  of  correct  discrimination  (PdO-  Pdj 
represents  the  fraction  of  1000  trials  for  which  the  algorithm  discriminates  correctly 
between  two  transients.  Pdi  curves  were  generated  for  SNR  values  of  10  dB  to  20  dB  for 
both  the  subspace  and  bispectrum  discriminators.  The  probability  of  incorrect 
discrimination,  Pfa,  was  computed  in  a  similar  manner  and  also  used  as  a  performance 
measure. 

Figure  22  shows  the  results  using  the  subspace  discriminator  where  the 
exponentially  decaying  sinusoid  is  present  at  sonabuoy  1  and  various  transients  are 
present  at  sonabuoy  2.  For  these  curves,  the  vertical  axis  is  labeled  "Probability"  because 
the  Pfa  and  Pm  are  shown  on  the  same  plot.  In  Figure  22,  the  Pdj  curve  for  the  exponential 
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Figure  22.  PD,  and  Pfa  (using  subspace  methods)  between  exponential  decaying 
transient  and  exponentially  decaying  transient,  CW  Pulse,  noise,  LFM  pulse  and 
whale  transient. 

decaying  transient  arriving  at  the  first  and  second  sonabuoy  is  shown.  It  can  be  seen 
from  this  curve  that  the  algorithm  discriminates  correctly  with  a  Pdj  of  0.9  or  larger  for 
SNR  values  greater  than  17  dB.    Figure  22  also  shows  the  Pfa  for  the  other  transient 
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arrivals  at  the  second  sonabuoy  (i.e.,  CW  pulse,  LFM  pulse,  or  the  whale  transient)  and 
the  Pfa  when  only  noise  (no  transient)  is  present  at  the  second  sonabuoy.  These  curves 
show  that  the  Pfa  is  very  low.  These  are  desirable  results  since  they  indicate  that  the 
algorithm  can  differentiate  between  the  transients  with  a  high  probability  of  correct 
discrimination  and  a  low  probability  of  false  alarm.  Results  for  the  other  transients  are 
summarized  in  Table  4;  the  corresponding  curves  are  given  in  appendix  A.  It  is 
important  to  note  that  the  results  for  the  whale  transient  are  based  on  a  longer  data 
observation  interval  than  the  others. 


Transient 
Arriving 

B 
U 

o 

Y 

2 

CW  Pulse 

Exponential 

LFM  Pulse 

Whale 

Noise 

BUOYl 

CW  Pulse 

PDi  >0.35  for 
SNR>19dB 

Pfa  <  0.08 

Pfa  <  0.01 

Pfa=0 

Pfa  <  0.01 

Exponential 

Pfa  <  0.02 

PDl>0.9       for 
SNR>17dB 

Pfa  <  0.002 

Pfa=0 

Pfa  <  0.021 

LFM  Pulse 

P&  <  0.03 

Pfa  <0.08 

PDi>0.8        for 
SNR>18dB 

Pfa  <  0.1 

Pfa  <  0.04 

Whale 

Pf3<0.2 

Pfa<0.18 

Pfa<0.18 

PDl>0.8     for 
SNR>18dB 

Pfa  <  0.19 

Table  4. 


PDi  results  using  subspace  methods. 


For  the  bispectral  discriminator,  the  threshold  Qt  (see  Eq.  3.44)  was  adjusted  by 
using  a  threshold  gain  to  improve  the  results.  The  gain  was  applied  to  Eq.  3.44  in  the 
following  way: 


Qr  = 


r       (N-\  - 

I>M-e-& 

\,=0 


N-2 


(5.  13) 


where  G  is  the  threshold  gain.     Threshold  gains  of  1,  2,  4  and  6  were  used. 

Figure  23  shows  the  discrimination  results  for  the  exponentially  decaying 
transient  using  the  bispectral  discriminator.  Again,  the  vertical  axis  is  labeled 
"Probability"  because  both  the  Pfa  and  Pdi  are  shown  on  the  same  plot.  For  Figure  23,  the 
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transient  arriving  at  sonabuoy  1  is  the  exponentially  decaying  sinusoid.  The  Pdj  for  this 
transient  is  shown  in  Figure  23.  The  other  curves  of  Figure  23  show  the  results  of  the  Pfa 
for  the  CW  pulse,  noise,  LFM  pulse  and  whale  transient  respectively.  The  symbols 
shown  in  Figure  23  are  tied  to  the  transient  type  and  are  consistant  with  what  is  used  in 
Appendix  B.  The  subplots  in  Figure  23  show  the  results  when  different  threshold  gains 
are  used. 
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Figure  23.        Pdi  (using  high-order  methods)  between  exponentially  decaying 

sinusoid  transient,  arriving  at  sonabuoy  1  and  the  CW  pulse,  exponentially  decaying 

Sinusoid,  noise,  LFM  pulse  and  whale  transients  arriving  at  sonabuoy  2.  Threshold 

gain  values  of  (a)  1,  (b)  2,  (c)  4,  (d)  6  were  used. 
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As  can  be  seen  in  Figure  23,  both  Pfa  and  the  PD,  decrease  with  increasing 
threshold  gain.  A  design  compromise  must  therefore  be  made  to  have  an  acceptably  high 
PDj  and  an  acceptably  low  Pfa.  For  this  thesis,  a  threshold  gain  of  4  was  chosen.  Table  5 
shows  the  discrimination  results  for  the  bispectral  discriminator  using  this  value  for  the 
threshold  gain. 


Transient 
Arriving 

B 
U 

o 

Y 

i 

CW  Pulse 

Exponential 

LFM  Pulse 

Whale 

Noise 

BUOY  1 

CW  Pulse 

PDi    >0.8    for 
SNR>17dB 

Pfa  <  0.2 

Pfa<0.01 

Pfa=0 

Pfa  <  0.001 

Exponential 

Pfa  <  0.3 

PDi   >0.98  for 
SNR>14dB 

Pfa<0.15 

Pfa  =  0 

Pfa  <  0.001 

LFM  Pulse 

Pfa  <  0.05 

Pfa  <  0.03 

PDi  >0.95  for 
SNR>16dB 

Pfa  =  0 

Pfa  <  0.001 

Whale 

Pfa=0 

Pfa=0 

Pfa  =  0 

PDj   >0.2  for 
SNR>20dB 

Pfa=0 

Table  5. 


Pdi  using  Bispectrum  using  a  threshold  gain  of  4. 


The  poor  discrimination  performance  of  the  whale  transient  given  in  Table  5  may 
be  attributed  to  the  relatively  short  observation  time  of  the  data.  Better  results  may  have 
been  obtained  if  longer  observation  times  were  used;  however,  it  was  not  practical  to  run 
a  large  number  of  simulations  for  these  longer  observation  times. 

C.         SUMMARY  OF  RESULTS 

The  TDOA  and  discriminator  results  are  summarized  in  Table  6  including  the 
advantages  and  disadvantages  of  both  algorithms.  As  has  been  previously  shown  in 
Figure  1 8,  the  bispectral  linear  phase  detector  tends  to  give  a  larger  number  of  correct 
TDOA  estimates  at  low  SNR.  On  the  other  hand,  the  subspace  algorithm  gives 
consistently  correct  results  (Pt  =  1.0)  at  high  SNR,  which  is  never  achieved  by  the 
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bispectrum  algorithm  for  SNR  in  the  range  of  8  dB  to  20  dB.     Therefore  the  more 
desirable  result  may  depend  on  the  application  where  the  algorithms  will  be  used. 

For  discrimination,  the  bispectral  method  gave  the  highest  PDi  results  when  a 
threshold  gain  of  4  was  used.  However,  the  Pfa  increased  with  increasing  SNR  (see 
Figure  23).  This  is  an  undesirable  result  since  even  transients  with  high  SNR  will  not  be 
able  to  be  separated.  On  the  other  hand,  the  subspace  discriminator  gave  lower  Pdj  values 
but  had  the  advantage  that  the  Pfa  remained  constant  for  all  SNR  values  tested. 


Algorithm 

TDOA 

Discrimination 

Bispectrum 

Best  Result:   SNR  >  lldB 
gives       a       Py       of      0.9 
(exponential  transient) 
Disadvantage:  A  Pj  of  1  is 
never  achieved  at  SNRs  < 
20dB 

Advantage:    High  Pi's   are 
achieved  at  low  SNRs 

Best  Result:  SNR  >  14dB  gives  a 
Pdj  of  0.98  (exponential  transient). 
Disadvantages:  Results  depend  on 
choice    of   threshold    gain.       Pfa 
increases  with  increasing  SNR. 
Advantages:   Can  achieve  a  high 
Pdj   at  relatively  low  SNRs.    Has  a 
low  Pfa  for  noise 

Subspace 

Best  Result:  SNR  >  17dB 

gives  a  PT  of  1 .  (exponential 

transient) 

Disadvantage:      High  Pt's 

are    only   achieved    at   high 

SNRs 

Advantage:     Pj's  of  1   are 

achieved. 

Best  Result:  SNR  >  17dB  gives  a 
Pdj  of  0.9  (exponential  transient). 
Disadvantages:      Computationally 
intensive.  Low  PDj  at  high  SNRs 
Advantages:     Low     Pfa     for     all 
transients 

Table  6. 


Summary  of  TDOA  and  discrimination  results. 
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VI.    CONCLUSIONS 

A.         THESIS  SUMMARY 

In  this  thesis,  we  have  developed  and  compared  two  algorithms,  namely  the 
bispectrum  and  subspace  linear  phase  detectors.  These  algorithms  were  developed  for 
the  purposes  of  transient  discrimination  and  TDOA  estimation,  in  order  to  be  used  as  part 
of  a  transient  tool  suite  to  aid  in  the  estimation  of  a  submarine's  position.  Chapters  III 
and  IY  of  this  thesis  provide  an  analysis  of  these  two  algorithms.  Two  performance 
measures  were  used  to  evaluate  the  algorithms,  namely  the  probability  of  correct  TDOA 
(Pj)  and  the  probability  of  correct  discrimination  (PdO- 

Chapter  V  discusses  the  simulations  and  transients  used  in  the  evaluation  of  the 
algorithms.  Of  the  four  transients  used,  it  was  found  that  the  whale  transient,  which  is  the 
longest  transient,  gave  the  worst  performance.  This  may  be  due  to  the  relatively  short 
data  observation  times  used  in  evaluating  the  algorithms.  It  was  found  that  longer 
observation  times  produced  better  results  for  this  transient;  however,  it  was  impractical  to 
run  a  large  number  of  cases  at  the  longer  observation  times. 

It  was  found  that  both  algorithms  had  advantages  and  disadvantages  as 
summarized  in  Table  6.  In  general,  it  can  be  said  that  for  TDOA  the  bispectral  linear 
phase  detector  gave  better  results  at  low  SNR  while  the  subspace  linear  phase  detector 
worked  better  at  the  higher  SNR. 

For  discrimination,  it  was  found  that  the  bispectral  discriminator  gave  higher  Pdj 
than  the  subspace  discriminator.  However,  the  Pfa  of  the  bispectral  discriminator 
increased  at  higher  SNR  while  the  subspace  discriminator  gave  a  constant  Pfa  for  all  SNR 
values  tested.  For  discrimination,  the  advantage  of  having  a  constant  Pfa  is  desirable. 
Therefore  the  subspace  discriminator  is  the  best  option  even  although  it  produced  lower 
Pdj  than  the  bispectral  discriminator.  As  discussed  in  Chapter  V,  there  are  design  trade 
offs  between  processing  speed  and  performance  that  need  to  be  made.  For  the  bispectral 
linear  phase  detector,  this  trade  off  is  in  terms  of  threshold  gain;  for  the  subspace  linear 
phase  detector,  this  trade  off  is  in  terms  of  correlation  matrix  size. 
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B.         FUTURE  WORK 

The  results  of  this  thesis  lead  to  some  interesting  topics  for  future  research.  These 
topics  include: 

(1)  The  effects  of  the  environment  on  a  transient.  In  the  evaluation  of  the  two 
algorithms  conducted  in  this  thesis,  no  work  was  done  to  determine  the  sensitivity  of  the 
algorithms  to  changes  in  transient  shape  due  to  the  environment.  To  do  this,  a 
propagation  model  must  be  used. 

(2)  The  effects  of  colored  noise.  In  the  simulations  additive  white  Gaussian  noise 
was  used.  This  is  a  poor  reflection  of  reality  since  ocean  noise  is  frequently  colored  [Ref 
22].  It  will  therefore  be  important  to  see  how  these  algorithms  perform  when  additive 
colored  noise  is  used. 

(3)  Operational  testing.  The  algorithms  were  tested  using  synthetic  data. 
Ultimately,  there  is  a  need  to  test  them  on  real  data  to  obtain  a  measure  of  their 
performance  and  applicability  in  the  field. 
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APPENDIX  A  PDI  AND  PFA  USING  SUBSPACE  METHOD 
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Figure  24.        PDi  for  the  CW  transient  arriving  at  the  first  sonabuoy. 


49 


0.9 


0.8 


0.7 


0.6 


ra    0.5 

.Q 
O 


0.4 


0.3 


0.2 


0.1  - 


t 


...*>v' 


0 


/• 


/ 


CW  Pulse,  P, 
fa 

Expo  Decay,  Pf 
Noise,  P, 


fa 


LFM  pulse,  P 


Di 


Whale,  P 


fa 


0^ 


_^_ 


_eb_ 


■&- 


-B- 
18 


12 


13 


14 


15 


16 

SNR(dB) 


17 


19 


20 


Figure  25.        PDi  for  the  LFM  transient  arriving  at  the  first  sonabuoy. 
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Figure  26.       PDi  for  the  whale  transient  arriving  at  the  first  sonabuoy. 
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APPENDIX  B 


PD,  AND  PFA  USING  BISPECTRUM  METHOD 
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Figure  27.       Pdi  (using  bispectral  methods)  between  CW  pulse  transient ,  arriving 
at  sonabuoy  1  and  the  CW  pulse,  exponentially  decaying  sinusoid,  noise,  LFM  pulse 
and  whale  transients  arriving  at  sonabuoy  2.  Threshold  gain  values  of  (a)  1,  (b)  2, 
(c)  4,  (d)  6  were  used. 
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Figure  28.        Pdi  (using  bispectral  methods)  between  LFM  pulse  transient ,  arriving 
at  sonabuoy  1  and  the  CW  pulse,  exponentially  decaying  sinusoid,  noise,  LFM  pulse 
and  whale  transients  arriving  at  sonabuoy  2.  Threshold  gain  values  of  (a)  1,  (b)  2, 
(c)  4,  (d)  6  were  used. 


54 


12      14     16 

SNR(dB) 
(a) 


12     14     16 

SNR(dB) 
(b) 


in 


—&, — ^v 

12      14     16     1*8   "  20 

SNR(dB) 
(c) 


1 

0.9 

[ 

- 

0.8  h 

- 

0.7 

- 

0.6 

1    0.5 

O 
°-    0.4 

- 

0.3 

- 

0.2 

- 

0.1 

- 

or 

*- 

0SE£iGSO£g    G-^-S"^ 

10 


12      14      16 

SNR(dB) 
(d) 


18 


20 


Figure  29.       PDj  (using  bispectral  methods)  between  whale  transient ,  arriving  at 
sonabuoy  1  and  the  CW  pulse,  exponentially  decaying  sinusoid,  noise,  LFM  pulse 
and  whale  transients  arriving  at  sonabuoy  2.  Threshold  gain  values  of  (a)  1,  (b)  2, 
(c)  4,  (d)  6  were  used. 
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